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ABSTRACT . 


This research develops computational methods for the digital 
simulation of continuous systems. These methods are applicable . to 
a wider class of systems and are more efficient than the methods used 
in currently available' simulation languages. Two fundamentally differ- 
ent approaches to simulation are distinguished and described, and the 
characteristics of the simulation which result from using these methods 
are investigated. 

Two serious limitations of available simulation languages are 
their inability to solve sets of implicit equations^'and to deal effi- 
ciently with stiff systems. Methods suitable for the efficient 
handling of stiff systems are developed within the framework of the 
two approaches to simulation. It is shown that these methods yield 
implicit equations for the overall simulation. In order to structure 
simulation languages to deal with implicit equations, sorting and 
iteration procedures for the efficient solution of these equations ■ are 
presented. 
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CHAPTER I 
INTRODUCTION 


1 . 1 Motivation 

As systems increase in size and complexity, simulation becomes 
an essential tool in their analysis and design. Simulation is used 
to corroborate theoretical derivations, to help in evaluation of alter- 
nate designs, to generate failure data without costly physical testing, 
to train human operators, to optimize design parameters, and for many 
other purposes. 

Simulation may be applied to either continuous -time or discrete- 
time systems. The variables of continuous -time systems (hereafter 
referred to as continuous systems) are capable of changing at any 
instant of time, whereas the variables of discrete-time systems (dis- 
crete systems) change only at discrete instants, and their values at 
all other points are of no interest. It is the simulation of contin- 
uous systems on the digital computer that will be considered in this 
thesis. A point worth noting, however, is that the modeling of con- 
tinuous systems for simulation on the digital computer necessarily 
involves the theory of discrete systems since the digital computer 
functions as a discrete-time machine. 

Prior to 1960, the simulation of continuous systems was carried 
out almost exclusively on the analog computer. Over the past decade, 
however, interest in the use of the digital computer has grown. 
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apparently beginning when the digital computer was first used to give 
check solutions for analog simulations. There are many reasons for 
the increased use of the digital computer for simulation, including 
greater accuracy, accessibility, no hardware setup, availability of 
mathematical functions and, most importantly, the existence of problem 
oriented simulation languages which allow the analyst to interface 
easily with the computer. Also, the growing availability of terminals 
and display devices is increasing the man-machine interaction, which 
was previously possible only with the analog computer. 

Simulation languages for the digital computer have grown to a 
high level of sophistication, and can simulate a large class of sys- 
tems. However, almost all languages are integrator oriented, that is, 
they separate integrators from the system and consider what remains as 
a functional block, which may be used as desired by a central integra- 
tion routine. While for many systems this is effective, if any sub- 
system does not lend itself to evaluation as desired by the central 
integration routine, then difficulties may arise. 

The goal of this research is to develop computational techniques 
both for the modeling of individual subsystems and for interconnection 
and coordination of the overall simulation. When incorporated into a 
general simulation language these techniques will extend the class of 
systems to which the language may be applied and increase its effi- 
ciency . 
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1.2 Historical Background 

The first digital analog simulator was reported by Selfridge in 
1955 [1.1], This was a digital computer program developed at the 
U.S. Naval Ordinance Testing Station, Inyokern, Calif., to accept a 
larger problem and provide more accuracy than was possible .on the 
available REAC (Reeves Electronic Analog Computer). Since that time, 
the field of simulation languages for digital computers has grown at 
a rapid rate. A history of simulation languages is given in several 
survey papers [1.2-1. 4], thus only that history which is necessary to 
motivate this research will be presented here. 

In their survey article, Clancy and Fineberg [1.3] state that the 
"essence" (in the metaphysical sense) of simulation languages may be 
parallelism, the apparent parallel operation of a serial digital com- 
puter. In many simulation languages, this apparent parallel operation 
is achieved by the analyst, who orders the input statements to the 
computer so that the proper computational sequence is followed at 
execution time. If one block follows another in the system under 
simulation, then the analyst must describe them to the computer in 
that order. Stein and Rose [1.5] in 1958 provided the theoretical and 
practical background needed to write a sorting routine, i.e., an' 
algorithm to process the input statements in whatever order they 
appear and to deduce the proper computational sequence from them. 
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Although overlooked by many later authors, a sorting algorithm is one 
of the keys to a useful language. A sorting algorithm is contained 
in most of the successful simulation languages, including MIDAS 
(Modified Integration Digital Analog Simulator) [1.6], MIMIC [1.7], 
PACTOLUS [1.8], DSL/90 [1.9] and most recently 360/CSMP [1.10]. 

Each of these languages, however, treats integrator outputs 
as known entities, and thus the sorting algorithm always begins the 
computational sequence at the integrator outputs. If the system being 
simulated has memoryless loops, or if any subsystem has been modeled 
by an implicit difference equation, then these sorting algorithms will 
fail. The existence of memoryless loops will be discovered by 360/CSMP 
and a message to that effect printed for the user, and the simulation 
will not be continued. The user may then break each memoryless loop 
with a special coding block, and, provided that there are no multiply 
imbedded loops, the program will attempt to solve the loop by a very 
simple iteration procedure. It is one purpose of the present investi- 
gation to develop a sorting algorithm which will accept memoryless 
loops and set up the necessary computational procedure to solve them 
by iteration. 

Even if the system being simulated contains no memory less loops, 
it is possible that they will be introduced through the modeling of the 
dynamic subsystems. In the literature, it has been shown [1.11] that 
numerical stability and efficiency are related to the use of implicit 
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numerical techniques. It will be seen that if subsystems are modeled 
using implicit techniques, memory less loops result. Thus, a strong 
motivation is provided for the development of techniques to deal 
effectively with the equations of memoryless loops. 

The numerical techniques referred to above lie in the domain of 
numerical analysis. Historically, the techniques of numerical analy- 
sis were the first to be applied to the discrete modeling of linear 
systems. After it became common to represent feedback control sys- 
tems by transfer functions, the Tustin [1.12] and other substitutional 
techniques for simulation [1.13-1.14] were developed. Other methods 
followed [1.15-1.17], Simultaneous with this work in simulation, many 
parallel concepts were developed in digital filtering [1.18]. Of 
particular interest in this respect is a paper by Steiglitz [1.19], 
in which an isomorphism between the discrete and continuous signal 
spaces leads to an identification of a discrete modeling method for 
continuous linear systems. In this thesis, the discrete modeling of 
systems will be approached from the point of view of numerical sta- 
bility and efficiency, and the relationship of the techniques devel- 
oped here to the work done in numerical analysis, digital filtering 
and simulation will be discussed. 
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1,3 Summary of the Thesis 

In Chapter II, a general model of a system as an interconnection 
of individual subsystems is given as well as a more precise definition 
of what is meant by simulation. Two approaches to the simulation of 
continuous systems are presented, and the modeling of linear and non- 
linear subsystems is discussed. Some limitations of conventional 
numerical methods when applied to stiff systems are shown. 

Dahlquist's [1.11] definition of A-stability is given, and the impli- 
cations of using A-stable methods are developed. 

The concept of A-stability is extended in Chapter III to the dis- 
crete modeling of linear systems, and a new technique for applying 
the bilinear transformation to irrational transfer functions is 
developed. 

A sorting algorithm for general simulation use is developed in 
Chapter IV, which sets up the computational sequence for solution of 
the problem by iteration. A presentation of iteration procedures 
for solution of nonlinear equations is also given. 



REFERENCES 


[1.1] R. C. Self ridge, "Coding a General Purpose Digital Computer 
to Operate as a Differential Analyzer," Proc. Western Joint 
Computer Conference (IRE) 1955 . pp. 82-84. 

[1.2] R. D. Brennan and R. N. Linebarger, "A Survey of Digital 

Simulation: Digital Analog Simulator Programs," Simulation 3 

no. 6, pp. 22-36 (December 1964). 

[1.3] J. J. Clancy and M. S. Fineberg, "Digital Simulation 
Languages: a Critique and a Guide," AFIPS Conference Proc. 27 
Part I, 1965, Spartan Books, Washington, D. C. , pp. 23-36. 

[1.4] J. C. Strauss, "Digital Simulation of Continuous Dynamic 

Systems: an Overview,” AFIPS Conference Proc. 33 FJCC 1968 , 

Spartan Books, Washington, D. C., pp. 339-344. 

[1.5] M. L. Stein and J. Rose, "Changing From Analog to Digital 
Programming by Digital Techniques," J. ACM 7 no. 1 (January 

i 960 ). 

[1.6] R. T. Harnett and F. J. Sanson, "MIDAS Programming Guide," 
Rept. No. SEG-TDR-64-1, Wright-Patterson AFB, Ohio (January 
1964). 

[1.7] F. J. Sanson and H. E. Petersen, "MIMIC-Digital Simulator 
Program," SESCA Internal Memo 65-12, Wright-Patterson AFB, 
Ohio (May, 1965). 

[1.8] R. D. Brennan and H. Sano, "PACTOLUS", AFIPS Conference Proc. 
FJCC 26, 1964, Spartan Books, Washington, D.C. 



8 . 


[1.9] W. M. Syn and R. L. Linebarger, "DSL/90-A Digital Simulation 
Program for- Continuous System Modeling," AFIPS Conference 
Proc. SJCC 1966 , Spartan Books, Washington, D.C. 

[1.10] System/ 360 Continuous System Modeling Program-User 7 s Manual 
(H20-0367) IBM Corporate Data Processing Division, 112 East 
Post Road, White Plains,- N. Y. 

[1.11] G. G. Dahlquist, M A Special Stability Problem for Linear 
Multistep Methods," BIT 3 no. 1, pp. 27-43 (1963). 

[1.12] A. Tustin, "A Method of Analyzing the Behavior of Linear 
Systems in Terms of Time Series,” J.IEE 94 pt. II-A, pp. ISO- 
142, (May, 1947). 

[1.13] A. Madwed,~ "Number Series Method for Solving Linear and Non- 
linear Differential Equations," Rept. No. 64445-T-26, Instru- 
mentation Laboratory, Massachusetts Institute of Technology, 
Cambridge, Mass., (April, 1950). 

[1.14] R. Boxer and S. Thaler," A Simplified Method for Solving 
Linear and Nonlinear Systems," Proc. IRE 44 , no. 1, pp. 89- 
101, (January, 1955). 

[1.15] M. E. Fowler, "A New Numerical Method for Simulation," 
Simulation 4, no. 5, pp. 324-330 (May, 1965). 

[1.16] J. S. Rosko, "Improved Techniques for Digital Modeling and 
Simulation of Nonlinear Systems," AFIPS Conf. Proc. SJCC 1968 , 


Spartan Books, Washington, D. C., pp. 473-481. 




9 . 


[1.17] 

[1.18] 
[1.19] 


A. P. Sage and S. L. Smith, . "Real-time Digital Simulation 
for Systems Control," Prop . IEEE 54, no. 12, pp. 1802-1812 
(December, 1966). 

F. F. Kuo and J. F. Kaiser, eds.. System Analysis by Digital 
Computer , John Wiley, New York, N. Y. , (1966), Chapter 7. 

K. Steiglitz, "The Equivalence of Digital and Analog Signal 
Processing," Information and Control 8, 455-467, (1965). 



CHAPTER U 


INTERCONNECTED SYSTEMS 


2,1 Introduction 

In this Chapter, two basic methods for the simulation of continuous 
systems are given. In one approach, the overall system modeling (OSM) 
method, the continuous input-output-state relations of the individual 
subsystems are combined to obtain one relation for the overall system. 
This combined relation is then digitally simulated. The limitations of 
this approach, when applied to subsystems of many different mathematical 
types, is shown. 

A second approach, the individual subsystem modeling (ISM) method, 
overcomes these limitations. The ISM method chooses a discrete model 
for each individual subsystem and interconnects these models to pro- 
duce the overall simulation. Techniques for modeling individual sub- 
systems are classified here. Most importantly, the stability and effi- 
ciency of the simulation for the different classes of modeling methods 
are discussed. Of particular interest is the use of A-stable tech- 
niques when dealing with stiff systems, i.e., systems with widely 
varied time constants. It should be noted that, for the OSM method, 
once the combined continuous system relation has been obtained', the 
completion of the method can be considered an application of the ISM 
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method to a system containing only one subsystem. Hence, the discussion 
of stability and efficiency for the ISM method is relevant to the OSM 
method as well. 


2.2 System Description 


In order to define exactly what is meant by a digital simulation, 
a mathematical model of a system as an interconnection of subsystems 
must be given. A collection of N subsystems is considered, each de- 
fined by an input-output-state relation 


y ± Ct) =A i [x i (t o )s v i» t3 

where y^(t) is the vector output of the i a subsystem, x^(t ) is the 

initial state and v^ is the vector input; is a function of the 

initial state x.(t ) and a functional of v. and t on the interval 

Ct ,t]. Only causal subsystems will be considered, i.e., for any two 
1 2 

input vectors v. and v. , if 

X X 

v^ = v? for t = t =.t., 
ii o 1 

then causality implies 


A i [x i (t o ) > v i ’ t] = i i [x i Ct o ) > v i > t] 



for all t Q , t^ and x^(t Q ) . 
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An interconnected system formed from these subsystems can be 
modeled as shown in Fig. 2.1. Letting 

v(t) = col[v 1 (t),v 2 (t), . .., v N (t)3 

and defining y( t) and x(t) similarly as the column vectors of outputs 
and states, an overall relationship for the interconnected system may 
be written as 


v(t) = M £(t) + P u(t) (2.1) 

where u(t) is a vector of external inputs. For simplicity, it will be 
assumed, with no loss of generality, that each component of a subsys- 
tem 1 s input vector is connected either to one component of another 
subsystem’s output vector or to one of the external inputs. In 
equation (2.1), this means that each row of the combined matrix 
CM P] has one and only one non-zero entry. The input-output-state 
relationships for all of the subsystems can be written collectively as 


Z(t) = A[x(t o ), v , t ] 


( 2 . 2 ) 


where 


A ^x(t ),v,t3 




WV’V 1 
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If equation (2.2) is substituted into equation (2.1), we have 

v(t) = MACx(t Q ) ,v,t] + Pu(t) (2.3) 

Since u(t) represents known external inputs, a solution of the system 

equations is the calculation of v in equation (2.3) and thus the 

determination of £(t) through equation (2.2). 

A digital simulation of a system of the form shown in Fig. 2.1 

is a numerical approximation to the solution of equation (2.3) at a 

n 

sequence of discrete time points t = t Q + ^ h. . The step size, 

j=l ^ 

h . , may be fixed or variable depending on the type of numerical tech- 
nique used. Approximations y n and/or v^ are desired such that 

Y n ~ £<t n ) 

- S< V 

2.3 The Overall System Modeling Method 

In the OSM method of digital simulation, the individual sub- 
system input-output-state relations are combined to obtain a single 
relation for the overall system. A discrete model is then chosen 
for this combined 'relation to produce the simulation outputs. If 
each subsystem is described by a set of normal form differential 
equations or by an equivalent block diagram, then each input-output- 
state relation is of the form 
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X.(t) = f i Cx i (t),v i (t),t] 

i = 1,2, ... N (2.4) 

y ± (t) = g i [x i (t),v i (t),t] 

If the output equation in (2.4) contains v^(t) as an argument, 
as shown, this represents an instantaneous feed-forward of information 
from input to output. If equations (2.4) are combined into a single 
vector equation, there results 


i(t) = F[x(t) ,v(t),t] 
£(t) = G[x(t),v(t),t] 


(2.5) 


where 



f 1 Cx 1 (t),v 1 (t),t] 


g 1 Cx 1 (t),v 1 (t),t] 

F[x(t),v(t),t] = 

f 2 [x 2 (t),v 2 (t),t] 

♦ 

and G[x(t),v(t) ,t] = 

g 2 [x 2 (t),v 2 ( t ),t] 


• 

• 




Substituting the combined output relationship of (2.5) into 
(2.1) yields 


v(t) = M G[x(t),v(t),t] + P u(t) 


( 2 . 6 ) 



16 . 


Equation (2.6) will generally be implicit in v(t) because of the 
instantaneous connections from inputs to outputs in the individual 
subsystems. Equations (2.5) and (2.6) are a canonical form for the 
overall system. Equation (2.5) represents N coupled differential 
equations and equation (2.6) are N auxiliary constraint relations. 

The implications of this canonical form for computer-aided analysis 
of nonlinear networks is discussed by Stern [2.1]. 

If equation (2.6) is implicit, then any of the large simulation 
programs will have difficulty. Implicit equations represent memory less 
loops in the original system, and of the more advanced simulation 
languages, only a few will attempt to solve these equations. If the 
loops are first found by the user, and if there are no loops imbedded 
within one another (equations must be uncoupled), then some of the 
programs will attempt solution by simple iteration procedure which 
converges only under strict constraints on the form of the functions 
involved. The sorting and iteration procedures developed in Chapter IV 
provide the necessary methods for solving equation (2.6) for very 
general systems. 

To continue with the OSM method, if the simulation is limited to 
subsystems with no direct connections from input to output, then 
equation (2.6) specializes to 


v(t) = M G[x(t),t] + P u(t) 
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and this may. be substituted into equations (2.5) to give 


*(t) = F[x(t) ,M G[x(t),t3 + P u(t),t] 
y.(t) = G[x(t),M G[x(t),t] + P u(t),t] . 


(2.7) 


This is the normal form representation for the overall system. Any 
standard integration techniques e . g. , Runge-Kutta, may now be applied 
to this combined set of equations. The available simulation languages 
[see Refs. [1.2-1J.0] of Chapter I] use essentially this approach, 
although equation (2.7) is never explicitly generated by the program. 
Instead, the program determines the proper computational sequence 
(sorting) so that evaluation of the subsystems in the prescribed order 
yields the necessary derivative values 1 

While any integration technique can be applied to equation (2.7), 
the available simulation languages are restricted to the use of ex- 
plicit methods, i.e. , methods which use only past values of input, 
output, and state to calculate the present output. The reason for 
this is that the sorting algorithms assume that integrator outputs are 
known quantities and thus are available for evaluation of derivatives. 
For implicit integration techniques, this is not true, because impli- 
cit equations must be solved at each time point to obtain the output 
or state. The limitations of explicit integration techniques, and 
hence of presently available simulation languages will be shown in the 
next section. 
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If ail the individual subsystems are described by normal form 
differential equations, then the OSM method can be used. If implicit 
equations of the form (2.6) are present, or if it is desired to use an 
implicit integration scheme, then techniques must be developed to solve 
the resultant implicit equations. The sorting and iteration procedures 
presented in Chapter IV are designed for this purpose. 

If, however, any subsystem is not in normal form, for example, if 
only its impulse response is given, then it may not be possible to use 
the OSM method, because a combined relation for the overall system may 
be difficult or impossible to obtain. In this case, another approach 
must be taken. When the subsystems are of different mathematical 
types, it is preferable to develop a discrete model for each of the 
individual subsystems and then to combine these models to obtain the 
overall simulation. This method, called the individual subsystem 
modeling (ISM) method, is described in the next section. 

2.4 The Individual Subsystem Modeling Method 

• The ISM method for digital simulation of continuous systems, as 
indicated in the previous section, is first to approximate each sub- 
system by a discrete model, and then to interconnect these individual 
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models in order to obtain a simulation of the overall system. In this 
section, techniques applicable to the modeling of general linear and 
nonlinear systems will be discussed. A detailed development of dis- 
crete models for linear systems will be left for Chapter III. 

The subsystem easiest to model is a memory less subsystem des- 
cribed by an input-output relation of the form 

y(t) = g(v(t),t) 

This equation is carried over into a discrete model of the form 


=' g <vV < 2 - 8 > 

If all subsystems are memoryless, e.g. , resistive electrical networks, 
then when the individual models are interconnected, the discrete ver- 
sion of equation (2.3) for the overall system becomes 


v = M G(v jt ) + P u 
— n — n n — n 


g, (v, ;t ) 
a l v l,n* n' 


6 <2n>V = 


s 2 (v 2,n ;t n> 


g N^ v N J n ;t n^ 


(2.9) 


where 
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A simulation requires the solution of the N coupled implicit equations 
in (2.9) at each time point, t^. This set of equations is usually 
sparse: therefore, in order to calculate the solution efficiently, 

advantage should be taken of the structure of the system [2.2]. The 
sorting and iteration procedures of Chapter IV are designed for this 
purpose. The modeling of memory subsystems, including all dynamic 
systems, is an area of wide interest and varied approach. It is 
possible, however, to form some classifications of the available 
techniques which help in determining the nature of the interconnection 
process once the individual models have been chosen. 

The general form for a discrete model of the input-output-state 
relation of a subsystem can be written either as 


y n “ p(y o> - ' * ’ y n ; V * * * ’ W ■ * • ’V 


( 2 . 10 ) 


or as 


y =r(x o ,...,x n ; V ...,Vt 0 V 


(2.11a) 

(2.11b) 


where, as in Section 2.2, y^ is the output vector, 
vector and v. the input vector at time t . 


the state 


Normally only a few of the arguments shown in the equations are 
actually present. The classification of the discrete models is deter- 
mined by the presence or absence of certain arguments as set forth in 
the following definitions. 
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Definition 2.1 : A discrete model of the form (2.10) which con- 

tains as an argument, or a discrete model of the form 
(2.11a) which contains x n as an argument, is called internally implicit . 

The presence of as an argument in (2.10) or of x n in (2.11a) 
indicates that the right hand side cannot be explicitly evaluated at 
time t . The word ” internally" is used in the definition to distin- 
guish an equation of the form (2.10) or (2.11a), which must be solved 
in order to update an individual subsystem, from other implicit rela- 
tions which may arise when all the individual models are interconnected. 
Subsystems of the following type generate these equations for the 
overall system. 

Definition 2.2: A discrete model is called externally memoryless 

if it is of the form (2.10) and contains v^ as an argument or if it is 
of the form (2.11) and either: 

i) contains v as an argument in (2.11b) or 

ii) contains v n and x n as arguments in (2.11a) and (2.11b) 
respectively. 

This terminology results from the following point of view: If 

the subsystem is considered to store all necessary past values of input, 
output, or state, then the appearance of v in equation (2.10) or 
(2.11) indicates that the output at time t^ is directly dependent on 
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the input at time t , which is similar to the discrete model of a 

n 7 

memoryless subsystem, equation (2.8). 

To illustrate these definitions, consider -the following example: 

Example 2.1: A set of normal form differential equations of the 

form 

x = f(x,v,t) 
y = g(x,t) 

can be approximated by many techniques. Euler integration. 


x = x , + hf(x , ,v . ,t t) 
n n-1 n-l ? n-1* n-1 


y - g(x ,t ) 

y n . n’ n 


is neither internally implicit nor externally memory less. Backward 
Euler integration. 


x = x , + hf(x ,v ,t ) 
n n-1 ^ n’ n J n 


y = g(x ,t ) 
J n n’ n 7 


however is both: internally implicit because x^ appears on both sides 

of the first equation, and externally memoryless by Definition 2.2ii. 
The second order Runge-Kutta formula 
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n 


n 


n 


= Vi + “(Vi’Vi-Vi 1 
= Vl + Nvi’Vl’W + ^VVV 3 

= g <vV 


is also externally memoryless, but is not internally implicit. 

An example of the fourth possibility, internally implicit but not ex- 
ternally memoryless, is not easily given for this example, since for 
most standard integration methods, the appearance of x^ is usually 
linked with the appearance of v . For general subsystems, however, the 
two definitions are independent. 

At the beginning of this section, it was stated that each sub- 
system would be approximated by a discrete model, and the individual 
models then interconnected to produce the overall simulation. In the 
original physical system, the connected variables are equal for all 
time. A discrete model of a continuous subsystem, as we have defined 
it in this section, produces output values, y , only at the discrete 
time points of the simulation, t , and values at any other time points 
are not available. Accordingly, when an output of one subsystem is 
connected to the input of another, equality is defined only at the 
points t . 
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To see how this discretization may affect the choice of numerical 
techniques used, consider a 4th order Runge-Kutta integration method 
for a set of normal form equations of the type considered in Example 
2.1. This method requires evaluation of the quantity 

M ( x n-i + ~F > v n-l/2 > ^ - I) 

where 

Vl/2 = V ( t n - I) 

and has been previously evaluated. If the subsystem were to be 
used alone, with v(t) an external continuous input, then v(t - h/2) 
would be available exactly. If, however, the input to this subsystem 
is the output of another subsystem, then only values v n _-^ and v^ are 
available, while the intermediate value v n _^/2 nat * Of course 
interpolation may be used to approximate the desired value, but in 
general this causes a loss in accuracy which negates the purpose of 
using a high-order integration technique. This discussion leads to 
the following definition. 

Definition 2.5: A discrete model is admissible for the ISM method 

if it is of the form (2.10) or (2.11), i.e., it only requires inputs at 
the discrete time points of the simulation. 
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In all that follows, only admissible models for the ISM method 
will be considered. 

In Example 2.1, it was seen that the choice of the numerical tech- 
nique to be used determines whether or not a subsystem is modeled as 
externally memoryless. The effect upon the overall simulation of 
modeling all memory subsystems as externally memoryless will now be 
determined. ' - 

The notation of equation (2.10) or (2.11) for an externally 
memoryless subsystem can be simplified if all arguments for past time' 
points are suppressed. It will also be assumed for the moment that the 
subsystem is not internally implicit. With these conventions understood, 
the input -output -state equation of an externally memoryless subsystem 
can be written exactly as that of an ordinary memoryless subsystem 

y = r(v ,t ) . 

•’n ^ n s n 

Equation (2.1) for the overall simulation then becomes 


v = M R (v ,t ) + P u 
— n — n 7 n — n 


( 2 . 12 ) 


where all variables are defined as in (2.7). Again, as in the case of 
a system composed only of memoryless elements, a sparse system of 
possibly nonlinear implicit equations must be solved. Had the 
dynamic subsystems not been modeled as externally memoryless, then this 
would have been the case only if there were loops of ordinary memoryless 


elements . 
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me special case where one or more subsystems are internally im- 
plicit will now be considered. An internally implicit subsystem has, 
by definition, an implicit equation associated with it. In order to 
update the subsystem, this equation must be solved. Two very different 
computational approaches to the solution of this equation are possible. 
The first approach assumes that this equation is solved internal to the 
subsystem, by its own solution technique, e.g. , iteration . From an 
input output point of view, the implicit equation does not exist, so 
that the previous discussion of externally memory less subsystems still 
applies, and equations of the form (2.12) result for the overall simu- 
lation. This approach ignores the interconnection of the subsystem with 
other subsystems as far as its internal implicit equation is concerned. 

A more efficient approach, however, is to consider all the implicit 
equations of the entire simulation simultaneously. These include both 
the equations of internally implicit subsystems and the equations which 
result when externally memoryless subsystems are interconnected. A de- 
tailed description of exactly how this is accomplished is given in 
Chapter IV. 

The remainder of this section will show why it may be desirable 
to model dynamic subsystems as externally memory less even though the 
added complexity of implicit equations results. The discussion first 
will be limited to subsystems described in normal form. It will be 
shown that for these subsystems the motivation for externally 
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memoryless modeling stems from a consideration of the stability and 
efficiency of numerical integration methods. For simplicity of 
notation, the system equations will be assumed to be in the form 

* = f(x,t), x(0) = x Q . (2.13) 

The input to -the subsystem, v, has not been shown, since, when model- 
ing an individual subsystem, the input is considered a known function 
of time. The output equation is also not shown, since it does not 
enter into the choice or use of an integration method for the state 
equation. The eigenvalues of the Jacobian, J = , are denoted 

X *2 ? • ♦ • 5 (arranged in order of increasing magnitude) and the 
linearized system is assumed to be asymptotically stable, i.e., 

Re(X ± ) < 0, i = 1, . . . ,N. 

In many applications it often happens that equation (2.13) is 
stiff, i.e., the stiffness ratio, p A JX N |/|\.J » 1. For these 
systems, the solution is smooth after a transient period during which 
there is a rapid variation. Examples abound in circuit theory, chemical 
kinetics, nuclear reactor calculation, and in many other fields. The 
majority of classical integration methods are unsatisfactory for stiff 
systems because the fast component of the solution continues to affect 
the numerical solution long after its effect in the' physical system 


has died out . 
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For most single-step methods, e.g., Runge-Kutta, the step size is 
limited by the largest eigenvalue of the Jacobian, J-, i.e., 

m * x |h X. | < d (2.14) 

1 1 jL 1 

where d' is a constant that is characteristic for the method. If 
this condition is not satisfied, then the solution will in general be 
numerically unstable, and hence useless. 

It is reasonable to assume that the solution interval, t, is 
proportional (by some factor c such as 3 or 10) to the inverse of 
the magnitude of the smallest eigenvalue, 

T = c[ m i n iX i j] 1 . (2.15) 

A lower bound on the number of steps necessary to compute the solution 
over the interval T is T/h where h satisfies (2.14), and thus from 
(2.15), T/h ~ p. If the spread of eigenvalues is large, then the num- 
ber of steps required will be excessive. For example, this bound can 
easily become arbitrarily large if some parameter, say a capacitor 
value in a linear R-C electrical network, becomes arbitrarily small. 
Methods having this limitation are called step length limited. 

Explicit multi-step integration methods have much the same problem 
when used to integrate stiff systems of differential equations. 
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The following example shows these difficulties for some particular 
integration methods. 

Example 2.2*: Consider the linear system of differential equations, 

* = Ax, x(0) = x Q 

with distinct real eigenvalues , 0 > X 1 > > . . . > . 

Since the analytic solution is 

x(nh) = ( e ^) x Q = x((n-l)h) 


and all eigenvalues are negative, the least that should be required is 
that 


For the explicit Euler, Heun or usual 4th order Runge-Kutta schemes, 
the approximate solution is 


x n = [M(hA)] n x Q 


(2.16) 


where 


1+q 

M( q) = 1+q+q 2 /2 l 

1+q+q 2 / 2 I +q 3 /3 \ +q^/4I 


Euler 

Heun 

Runge-Kutta 


*This example is taken from a short paper, "Stiff Systems of Differential 
Equations," by J. S. Rosenbaum, prepared for a course in nonlinear 
systems at Columbia University, (Fall, 1968) 
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Since the eigenvalues of A are distinct, there exists a matrix P 
such that 

P" 1 AP = diagC^, . . . , X r ] = A . 

Therefore, the difference equation (2.16) becomes 

x = P~ 1 [M(hA)] n Px 
n o 

which implies that 

lim x n = o iff |M(hX^)| < l i = 1, , N . 

n - oo 

In particular, for Euler's method, there results 

x = x . + h f(x t , t) 
n n-1 v n-1* 

or x = P _1 (l + h A) P x 

n v n-1 

therefore, stability requires that 

|1 + h X | <1 i = 1, . . . , N 
or |h X.| < 2 


For Heun' s method the stability condition is that |h X^| <2 and for 
the Runge-Kutta method the stability condition is jh X^j <2.78. 
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To demonstrate the difficulties that these conditions imply for 
stiff systems, consider the following second order example, 

v = Dv + C , v(o) = v Q 


where 


D = 

-1 

0 

and 

C = 

2 


0 

- 1 _ 



_2_ 


For this equation, v^(t) = v^( t) = 2(1- e and the stability conditions 

for Euler (h < 2), Heun (h-< 2) and Runge-Kutta (h < 2.78) place reason- 
able limits. on h. However, consider the system 


w = A w + C 


w(o) = 


-.1 

0 


where 



-500.5 

499.5 


2 

A = 



and C = 



499.5 

-500. 5_ 


_ 2_ 


for which the solution is 


w x = 2(1 - e -t ) 
w 2 = 2(1 - e" t ) 


, -lOOOt 
.1 e 

.1 e 


-lOOOt 
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For t > .02, .1 e “ 1000t <2.5 x 10~ 10 , or v « w, However, the eigen- . 
values of A are -1 and -1000, which requires that h be 1000 times 
smaller for the w equation than for the v equation. 

In the previous discussion, it has been seen that the efficiency 
of numerical integration methods when applied to stiff systems is 
closely related to the stability of the numerical solution. ' In many 
applications, it is stability rather than accuracy considerations which 
necessitate the use of a very small step size. Dahlquist [2.3] has 
given a definition of stability for numerical methods which leads to 
the identification of integration methods suitable for stiff systems. 

Definition 2.4: A numerical integration method is called A-stable 

if all solutions tend to zero as n -1 “ when the method is applied with 
fixed positive h to any differential equation of the form 


x = qx 

where q is a complex constant with negative real part. 

This definition states that for A-stable methods, the transient 
solution eventually goes to zero, and also that the transient solution 
does not affect the step size. 

It can easily be shown that the usual explicit 4th order Runge- 
Kutta method is not A-stable, for the Runge Kutta method yields a 
truncated Taylor series when applied to equation (2.15), i.e., it 
gives the sequence 
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[1 + qh + (qh) 2 /2l + (qh) 3 /3l + (qh) 4 /4 T .] n n = 0,1,2,... 


and this does not tend to zero everywhere in the left half plane, 

Re(q) < 0. Runge-Kutta methods are not admissible for the ISM method, but 
they can be used in the OSM method. If the OSM method is to be used 
efficiently for stiff systems, then, since explicit Runge-Kutta methods 
are not A-stable, implicit Runge-Kutta methods must be used- [2.4]. 

Implicit Runge-Kutta techniques give externally memoryless discrete 
models, because they always use the input at the present time. This , 
as shown previously, results in implicit equations of the form (2.12), 
which must be solved at each time point. 

A linear multi-step integration method has the form 


+ ax . = h(b. f 
o n-k s k n 


+ ... + b f .) 

o n~k' 


where a., b. are real constants, i = 0,1,2,... k, and f = f(x ,t ). 

If the vectors x q,x 1 , . . . and v^ , v^ , . . . , v^_^ are given, then 

x^,x^ + ^,... are computed recursively through (2.17). If the method 
is implicit, i.e., b^ ^ 0, then equation (2.13) must be solved for x n 
at each time point. In the notation of this section, b^. ^ 0 also indi- 
cates that the technique is externally memoryless. Since only values 
at the discrete points are required, linear multi-step methods are 
admissible for the ISM method, and hence the discussion to follow is 
relevant to both the ISM and OSM methods. 

Dahlquist [2.3] proves the following theorem, presented here 
using the terminology of this Chapter. 
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Theorem 2,1: An A-stable linear multi-step method must he inter- 
nally implicit and externally memory less, i.e., ^ 0. 

Thus, in order to achieve A-stability of the simulation using 
linear multi-step methods, externally memoryless models must be used. 
It will be seen in Chapter III that this same motivation for exter- 
nally memoryless modeling exists in the simulation of linear rational 
and irrational transfer functions. 

To carry the discussion of linear multi-step methods further, 
it is necessary to define the order, p, of an approximate method. 
Assume that the values x n ,X-,,...,x are exact, i.e., x. = x(t. ), and 
that the value x n+1 is calculated using the method in question. Ex- 
panding both and x(t in a Taylor series about the point 

t = t . there results 
n 7 

x (t n+ i) = x(t R ) + h /li X 1 (t n ) + h 2 /2I x’ ' (t n ) + ... 

and 

x n+1 = x(t n ) + h/li c 1 + h 2 /21e 2 + ... 

The highest exponent of h for which corresponding terms in the two 
series are equal is called the order, p, of the method, i.e., 

e l = X'CV,..., e p = x (P) (t n ), but c p+1 f x (p+1) (t n > 
Dahlquist proves the following theorem. 
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Theorem 2.2: The order, p, of an A-s table linear multi-step method 

cannot exceed two. The smallest truncation error is obtained with the 
trapezoidal rule, 

x n = Vl + h/ 2 (f n-l + f n > 

A-s table linear multi-step methods for the solution of stiff 
systems of differential equations have been discussed in the litera- 
ture [2. 5-2. 7] and Liniger and Willoughby [2.5] give several examples 
of the dramatic increase in efficiency that is possible using these 
methods . 

The previous discussion has shown that externally memoryless 
modeling is useful in order to achieve an efficient simulation of 
stiff systems. Another motivation for the use of externally memoryless 
models exists when the waveforms of the system being simulated have 
discontinuities or discontinuous derivatives. In this case, the use 
of the present value of the input will tend to force the solution to 
react more quickly to the effect of the discontinuity. For example, 
if the input to a system in normal form (2.13) is a ramp, beginning 
at = T > t Q , then, if an explicit method is used, the next value 
x n+l calcu - a1:e ^ using only past values of input, which are iden- 
tically zero. An externally memory less technique, however, uses the 
value of the input at time t n+ ^(nonzero) , and hence calculates x n+ -^ 
using the fact that the input has changed. 
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In the next section, an example is presented which shows the 
benefits of externally memoryless modeling for systems with discon- 
tinuous waveforms. This example consists of a transmission line with 
nonlinear loads and thus represents two coupled subsystems of differ- 
ent mathematical types. For this reason, it is necessary to use the 
ISM method . 
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2.5 Transmission Line Example 

In this section, an example is presented which illustrates the use 
of externally memoryless modeling to achieve accurate simulation for an 
interconnection of two systems of different mathematical type. The 
problem considered is shown in Fig. 2.2, and consists of a step current 
source shunted by a parallel R-C network feeding a lossless transmission 
line. The transmission line is terminated in -a nonlinear load. As 
suggested by Wing [2.8], we may model this network as two interconnected 
subsystems, one representing the source and load, and one representing 
the transmission line. Fig. 2.3 is a block diagram of this representa- 
tion. 

If the transmission line is regarded as a two-port, the terminal 
relation may be written as 

i(t) = h(t)*v(t) 

where i(t) = colCi^t), i 2 (t)] and v(t) = colCv-^t), v 2 (t)] , 
h(t) is the impulse response of the network and * denotes convolution. 
Thus, for this system, v(t) is regarded as the input and i(t) as the 
output. 

For the lossless transmission line, h(t) is found to be 

00 

' h ll = h 22 = s(t) + 2 I 8 (t-2k) 

k=l 

CO 

h 12 = h 21 = ” 'Z 6 (t-2k+l) 

k=l 
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FIG. 2.3 BLOCK DIAGRAM OF TRANSMISSION LINE PROBLEM 
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Since the impulse response contains only impulses, a numerical approxi- 
mation to the convolution of h(t) with v(t) becomes a simple weighted 

$ 

sampling of present and past values of v. The impulse at the origin 
in h . [ n and hg 2 samples the present value, v^, and hence this model is 
externally memoryless. If the transmission line were not lossless, a 
more complex h(t) would -result, with a corresponding increase in the 

t 

complexity of the discrete model. In general, however,- the model 
would still be externally memoryless. 

The discrete model of the lossless line can be reduced to 


where I ^ depends only upon the past values of the input, v. The 
sampling involved in calculating I^_^ at each step occurs because it 
takes two time units for a signal to travel down the line, be re- 
flected, and return. 

The lumped elements will be described by a set of normal form 
differential equations 


v = f(v,i) 

which for the network given is 


( 2 ‘. 18 ) 


f 

<• 

H 

l 


" (2 - v 1 /R 1 - i 1 )/C i 

_ V 2 _ 


(- g^) - i 2 )/C 2 


4 In Chapter III it will be shown that this is not necessarily the best way 
to approximate h(t), but it will be used here' for simplicity. 
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Here, i is now considered the input, and v the output. The trapezoidal 
rule will be used as a discrete model for the lumped elements, but 
modified, as suggested by Liniger and Willoughby [2.5], to include 
one Newton-Raphson iteration. 

The trapezoidal rule applied to equation (2.18) gives 


v. 


n 


= v , + ~ [f , + f ] 
n-i. i n-x n 


where f R = 

Letting z = v , this can be written 

Q(z) = z - - Vl = 0 

and taking one derivative with respect to z, 

Q J (z) = I - | f z (z,i n ) 
where f z = ^ . 

Taking one Newton-Raphson iteration using z q as the initial point gives 

z 1 = z° - [Q 1 (z)]" 1 


. Q(z°) . 



42 . 


Taking only this first iteration, letting z q = v^_^, there results the 
explicit integration formula 


v = v , + [I - 
n n-1 


h 

2 


f ( v 
v^ n- 


l ,a, n )] 


-1 


2*" f n-l + ^n-l’ 1 ^-* * ^ 2 * 19 ^ 


This formula is the discrete model for the lumped elements. It should 
be noted that the model is externally memoryless due to the presence 
of the input value i , but it is no longer internally implicit because 
the Newton -Raphson iteration has removed v as an argument on the right 
hand side. To simplify notation, equation (2.19) is written symbolic- 
ally as 


v 

n 


F(i. 


v 


n’ n-1 5 n-l 


( 2 . 20 ) 


Since both subsystems are modeled by externally memoryless 
techniques, it is expected that implicit relationships will exist when 
the individual models are interconnected. The interconnection has in 
essence been accomplished already, through the use of identical variable 
names, i and v, when modeling each subsystem. Therefore, the general 
equation for interconnected models, equation (2.1), becomes for this 
example 




~ 1 V 4- I ~ 

n n-1 

V 

n 


F(i ,i -j , v t ) 

L n’ n-1* n-1 _ 


( 2 . 21 ) 
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Solution techniques for the general case of equation (2.1) will be 
left for Chapter IV, but so as to be able to continue with the example, 
a simple iteration procedure will be used here. If the second equation 
in (2.21) is substituted into the first, there results 


i 

n 


F(i. 


n 


,i _ ,v ,) 
9 n-1 9 n-1 


+ I. 


n-1 


( 2 . 22 ) 


Since i , , v , , and I , are known quantities at time t , this 
n-1* n-1’ n-1 H n 5 

equation is of the form 


x = 9(x) . 


An obvious iterative solution technique is 


k+1 . k. 

x = cp(x ) 


which will be shown in Chapter IV to converge if 


H II ^ C < 1 


where ||*|| denotes the matrix norm ||A|[ A max Z |a..| . 

i 3 ^ 

This method is used to solve equation (2.22), resulting in 


i k+1 = F(i k , i . , v n ) + I i . 
n v n 9 n-1 9 n-1 n-1 
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At each time point t , an initial estimate i° must be chosen. The use 

n’ n 

of the previous value, i° = i , was first tried, but it was found that 

’ n n-1 5 

if a linear prediction was made, the total number of iterations neces- 
sary for convergence was reduced. Hence, 


i° = 2i .. - i _ 
n n-1 n-2 


Each time a reflected wave reaches one end of the line, discon- 
tinuities in the slopes of the variables occur. The use of externally 
memory less modeling for both the transmission line and the lumped 
elements helps in handing these waveforms, as explained in the previ- 
ous section. During calculation, it is seen that the iteration count 
in solving equation (2.22) is higher for the time points where the 
discontinuities occur than for time points where the waveforms are 
well behaved. 

The transmission line example discussed above was run on a time- 
shared SDS/940 computer for both linear and nonlinear loads. For the 
first simulation, R^ = IK, = lpf, R 2 = 2K, and C 2 = Ipf. Fig. 2.4 
is a plot of v- L ,v 2 »i 1 an ^ a 2 ^ or a ste P s i ze °f 1/32 ns. The iteration 
at each step was stopped if the difference between the predicted and 
calculated values was less than an iteration tolerance, £. 

Since with linear loads the entire problem is linear, it is 
possible to compute the exact values using Laplace transform techniques 
and working with the incident and reflected waves. These values agree 
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very closely with the computed results, the error being too small to 
show in Fig. 2.4. 

To determine the effect upon solution accuracy of choosing differ- 
ent values for the iteration tolerance, 6, and the step-size, h, 
several runs for different values were made. These results are 

summarized in Table 2.1. It is seen that reducing € beyond a certain 

* 

value simply increases the iteration count with no increase in accuracy. 
As would be expected, this value of € is approximately the same as the 
truncation errors introduced in the numerical approximation techniques. 

The linear resistor Rg is now replaced by a diode modeled by the 
diode equation 

v/v 

i = I 0 (e - 1) 

with I Q = .005 ma and v T = .026 v. A plot of the results for h =1/64 
ns appears in Fig. 2.5 and Fig. 2.6. It is not possible to obtain 
analytic results, but the computed results approach the correct 
steady-state values, and several runs with smaller step-sizes indicate 
that these results are approximately of the same accuracy as achieved 
for the linear example. 




I.Or- v, .VOLTS 



2 * 


VOLTS 
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TABLE 2.1 


Statistics For Transmission Line Example With Linear Loads, 


Step 

Size 

Iteration 

Maximum 

Error 

Average 

Error 

Iterations/ 

(ns) 

Tolerance 

(X 10 ~ 3 ) 

(X 10 -3 ) 

Step 

1/32 

• 

o 

H 

4.7 

1.51 * 

1.26 

1/32 

.005 

2.1 

.88 

1.45 

1/32 

.001 

1.7 

.32 

1.89 

1/32 

.0005 

‘ 1.6 

.34 

2.01 

1/32 

.0001 

1.4 

.32 

2.34 

1/32 

K set ss 2 

1.5 

.33 

2.00 

1/64 

.005 

2.1 

.599 

1.10 

\ 

1/64 

.001 

.47 

.172 

1.51 

1/64 

.0005 

.42 

.091 

1.72 

1/64 

.0001 

.38 

.081 

2.00 

1/64 

5 X 10“ 5 

.37 

.081 

2.07 

1/64 

1 X 10’ 5 

.35 

.076 

2.33 

1/64 

K set = 2 

.37 

.079 

2.00 

1/128 

.0005 

.260 

.077 

1.24 

1/128 

.001 

.099 

.021 

1.77 

1/128 

5 X 10" 5 • 

.094 

.018 

1.89 

1/128 

1 X io" 5 

.092 

.019 

2.04 

1/128 

5 X 10" 6 

.090 

.020 

2.10 

1/128 

K set ss 2 

."090 

.019 

2.00 
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CHAPTER III 

MODELING OF LINEAR SYSTEMS 


5 . 1 Introduction 

In this Chapter, general methods of A- stable modeling of linear 
fixed systems are presented. It is shown that A-stable techniques 
result in externally memoryless models, and hence the question of 
solution of implicit equations for the overall simulation again 
arises . 

A-stable techniques are applied to the irrational transfer func- 
tions that are typical of distributed systems. An example is given in 
order to illustrate the difficulties which may be encountered when the 
more conventional approaches to modeling are used. 
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3.2 A-stable Modeling of Fixed Linear Systems 

In Chapter II , a general formulation of a system as an intercon- 
nection of subsystems was given, and some general approaches to 
modeling of the subsystems were discussed. If each subsystem is now 
considered to be linear and fixed (but not necessarily lumped), then 
the following problem can be considered. Let the original intercon- 
nected system be stable. It is possible that the original system may 
contain some subsystems which are unstable, even though the overall 
system is stable. The problem is to determine a class of discrete 
modeling methods, which, if used to model each individual subsystem, 
gives A-stability for the overall simulation when the individual models 
are interconnected. Application of A-stable modeling to the unstable 
subsystems would generally produce unstable models. . What is required 
here, however, is that the interconnection of all the models, both 
stable and unstable, produce an A-stable overall model when the ori- 
ginal overall system is stable. It will also be required that the ' 
model be consistent, i.e. , a reasonable approximation of the original 
system. This problem of A-stable modeling will now be more precisely 
defined using some basic definitions and theorems (Theorems 3.1 - 3.4) 
from the theory of linear systems (see e.g., [3.1-2]. 

Definition 3.1 : A linear fixed causal continuous system, S, with 

input vector u and output vector y is stable if its zero-state 
response remains bounded for any bounded input. 



If F(s) represents the transfer function matrix of S, with input and 
output defined as above, and if f(t) is the impulse response of S , 
i.e., f(t) is the inverse Laplace transform of F(s), then the follow- 
ing well known theorem applies. 

Theorem 5.1 : S is stable if and only if each component of the im- 

pulse response, f. .(t) satisfies the condition 

J 

00 

J |f ij .(t)|dt < ~ . 

o 

It will prove useful later to have a necessary condition for stability 
in the frequency domain. 

Theorem 3.2 : A necessary condition for the stability of S is that 

i 

each component of the transfer function matrix, F^_.(s) , be finite for 
Re(s) 2= 0. 


Proof : 




-T,.. 

J. n 




(t)e st dt 


From this, there follows 


03 

|F i;j (s)| * J |f ±j (t)|e 


-at 


R e (s) * 0 


dt 
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“(Jt 

where c = Re(s). If Re(s) ^ 0, then e § 1 , and therefore 

CO 

|Fi .(s) 1 jf (t)|dt Re(s) ^ 0 . 

J o 3 

The right hand side is finite by Theorem 3.1, which completes the 
proof. Note that this theorem applies to irrational as well as 
rational transfer functions. 

Since an identification must be made between stability of a con- 
tinuous system and stability of its discrete model, the following 
definition and theorems for stability of discrete systems are given. 

Definition 3.2 : A linear fixed causal discrete system, D, with input 

vector u and output vector y, is stable if its zero-state response 
remains bounded for any bounded input sequence. 

If G(z) represents the transfer function matrix of D, and if g(n) is 
the response to the input sequence 1,0, ..., i.e., g(n) is the inverse 
z-transform of G(z), then the following theorem applies. 

Theorem 3.3 : D is stable if and only if each component of the zero- 

state response, g^(n), satisfies the condition 

00 

Y lsr ± j<n) | < 00 . 

n=o 

It is interesting to note that Theorem 3.1 does not require | f . • C t ) | 0 

"*■ J 

as t w but Theorem 3.3 does require Jg^(n)| ■* 0 as n 
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A necessary condition for stability in the frequency domain will 
be useful later. 

Theorem 5.4 : A necessary condition for the stability of D is that 

each component of the transfer function matrix, G..(z), be finite 
for | z| ^ 1. 

The proof follows that of Theorem 3.2 and makes use of Theorem 
3.3, and so will not be given here . 

A class of discrete models for continuous systems will be defined 
by a mapping from the s -plane to the z-plane. For each subsystem 
transfer function F^(s) let 


Fj(z,h) = F i (g(z,h)) 

where the mapping s = g(z,h) is chosen so that the resulting discrete 
model F^(z,h) has the following properties: 

(1) Consistency 
Let 

f i (n ’ h) = 2ij I FjC^tOz 1 *- 1 dz 
C 

where C is a circle enclosing all the singularities of F^(z,h)z n_1 , 

i.e., f^(n,h) is the inverse z-transform of F^(z,h). The discrete 
.model is then said to be consistent if 
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lim 

n-* 40 

h-*o 

nh=t 


h 


f i<« 


except possibly at discontinuities. 
Rewriting, it is required that 


lim 

n-*° 


n 


! i( n > k) 


= qct) 


(2) A-Stability 

If F^(s) is a stable transfer function, then F^(z,h) must be 
stable for all h ^ 0, i.e., the discrete model is stable for all 
positive values of the step-size. 

The second property implies the following theorem. 

Theorem 3.5 : If s = g(z,h) is an A-stable mapping, i.e., it satisfies 

condition (2) above, then 

Re(s) < 0 implies |z| < 1 . 

Proof ; F^(s) stable implies Re(s) <0 if F^(s) is not finite, from 
Theorem 3.2. F^(z,h) stable implies |z| < 1 if F|(z,h) is not 
finite from Theorem 3.4. Hence if s = g(z,h) satisfies condition (2), 
A-stability, then 


Re(s) < 0 implies |z| <1 , 
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Thus it is seen that the left half s -plane must map into the interior 
of the unit circle in the z -plane. 

It is obvious that if the same consistent A-stable transformation 
is made for each subsystem, then the overall model, formed from the 
interconnection of the individual models, will also be consistent and 
A-stable. It will be shown later ‘ however , that applying different 
transformations to subsystems, each of which is consistent and A- 
stable, will not lead in general to a consistent A-stable overall 
model. 

It now remains to determine functions g(z,h) which satisfy both 
conditions (1) and (2) above, that is, they are both consistent and 
A-stable. If the search is. limited to rational functions of z, then 
it is possible to derive certain characteristics of the allowable 
transformations. Consider a k-th order transformation of the form 


g( z jh) 


1 a(z) _ 1 
h * b(z) ~ h 



+ a o 
+ b Q 


(3.1) 


where a(z) and b(z) have no common zeros. If the system being 
modeled is of the form 1/s, then a transformation (3.1) is equivalent 
to using a linear multi-step integration formula. Dahlquist [3.3] 
has worked with A-stable integration formulas of this type, and it is 
possible to follow a similar development here. The following theorem 
establishes necessary and sufficient conditions for a transformation 
of the form (3.1) to be A-stable. 
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Theorem 3.-6 : A k-th order transformation s = g(z,h) of the form 

(3.1) is A-stable if and only if g(z) = a(z)/hb(z) is analytic for 

\z\ * l. 

Proof : If s = g(z,h) is A-stable, from Theorem 3.5, 


Re(s) < 0 implies jz| < 1 


or 


|z| ^ 1 implies Re(s) s 0 . 

For a transformation of the form (3.1), 


s = rtfi) if b < z) * 0 • 

If z 1 is a zero of b(z), then a(z^) ^ 0 since a(z) and b(z) have 
no common factors. In the neighborhood of z^, 

hb^T - ‘ c < M i> c * ° 

for some positive integer m. But Re[a(z)/hb(z) ] ^ 0 in a whole 
circle around z ^ if | ^ 1 for stability. Hence b(z) ^ 0 if 
|zj ^ 1, i.e., a(z)/hb(z) is analytic for jz| ^ 1 . 

Theorem 3.7 ; A k-th order transformation of the form (3.1) is not 
A-stable if b^ = 0. 
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Proof : For some integer m, m ^ 1, 

b(z) « c Z k ' m , when z -* “ , c ^ o . 
k 

However, a(z) a^z , ^ 0 . Hence 

rj d z m , where d ^ 0 , m ^ 1 . 

This however, contradicts the condition Re(s) s 0 for |zj ^ 1 

Theorem 3.8 has important implications for the modeling" of 
linear systems. For A-stable transformations, Theorem 3.8 implies 
that the resultant model is externally memory less, as shown below. 
Consider the initial value theorem for the z- transform 


fiCo^h) = 1 ™ F.[g(z,h)3 
For g(z,h) of the form (3.1), with b^, ^ 0 


lim 

Z -H» 


g(2,h) 


i 

h 



hence 

_ = F i[ K ' ^ ] 

and therefore f(0,h) will in general have a nonzero value. The use of 
the discrete model can be represented by the convolution 
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n 

y n = I f i«» h) Vj. 

j=o 

where y is the subsystem output and v is the subsystem input. 

The leading term in the convolution is f^(0,h) v^, and hence the 
present value of the input appears on the right hand side of the 
discrete input-output .relation. By Definition 2.2, the model is 
therefore externally memoryless, and implicit equations are expected 
for the overall interconnected model. It is important to note that 
these implicit equations result specifically from the A-stable 
modeling. 

As an example of a transformation of the Form (3.1), consider 
the well known bilinear transformation, 



z-1 

z+1 


This is both A-stable and consistent as is shown below: 


(3.2) 


(1) Consistency 

For the bilinear transformation applied to a transfer function 
F(s), consistency requires 


lim 

n-*° 


I %("> I) ■ f <« 


but 


lim n f ^ _ lim 1 X P z-ll n-1 , 

n— t f i< n > h) “ n -* 00 2ttJ J F ±ln * ITTJ 2 dz 
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where C is a circle centered at z=0, with radius r > 1 since the 
integrand has all singularities inside the unit circle. The integral 
along C in the z-plane can be transformed into an integral in the 
s -plane. 

If the inverse transformation 


z 



- (l - f^)(h/2) - (l + ^)(-h/2) 

dz = g ds 

(1 - sh/2; 


h ds 



is made, there results 


lim JL_ 
n-*° 2rrj 


t 

C ! 




where C T is a line parallel to the imaginary axis and an infinite 
semicircle. This can be reduced to 


and as n - * 30 


lk J/ 

n -«a 2 rrj 


1 


1 


st 

2n 


.-in 


st 

2n 


i 


F(s) 


1 + Sr 


st 

2n 


1 - 


st 

2n 


n 




5 


_i_st 

e . Thus, there results in the limit 
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1 

2TTj 


F(s) e st ds = f(t) 

v 

C 1 


and consistency has been proved. 


(2) A-Stability 


z = 


n sh 

1 + T 

, sh 
1 ~ 2 


1 ah . bh 

1 + — + 3 — 


ah 

2 


1 - — - J 


. bh 
2 


where s = a + jb. Clearly, if a < 0, |z| <1 for all h > 0. 


A system described by a rational transfer function can always 
be put into normal form, hence for all rational transfer functions, 
the following theorem of Dahlquist [3.1] applies. 


Theorem 3.8 : The order of an A-s table technique cannot exceed 2. The 

smallest truncation error is achieved with the bilinear transformation 
(trapezoidal rule). 

The bilinear transformation also has the desirable characteristic 
that rational transfer functions of order p transform into discrete 
models of the same order. 

The bilinear transformation has been used in many different areas 
and approached in many different ways. In numerical analysis, it 
represents the generating function of the trapezoidal rule, which was 
one of the first techniques to be used for the integration of 
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differential equations. In the field of digital filtering, it is used 
to obtain a discrete filter from a given continuous filter, with a 
compensation usually made because of a nonlinear warping of the fre- 
quency scale E3.4]. This nonlinear warping occurs because if 

s = £ fezi) 

s h \z+l / 

s^t 

and if z = e , then 


s 


= $• tanh 
h 


/ s i h 

\ 2 y 


\ 

t 


Working in digital signal processing, Steiglitz [3.5], has defined an 
isomorphism between the space I^C- 00 , 03 ) of continuous signals and the 
space of discrete signals. He then shows that this isomorphism can 
be used to associate a discrete filter with every continuous filter 
through the bilinear transformation. 

The bilinear transformation is not new to simulation, where it 
is called the Tustin [3.6] method. In the application of this method 
(to rational transfer functions), however, an arbitrary unit delay is 
inserted into every feedback loop to prevent the occurrence of implicit 
equations. The approach presented in this chapter does not insert any 
delays, for this in general destroys the A-stability of the model. The 
bilinear transformation will be considered further in the next section 
in order to explore the implications of A-stable modeling, especially 
when applied to irrational transfer functions. 
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3,3 The Bilinear Transformation 

The bilinear transformation is equivalent to the trapezoidal 
rule, which was discussed in Chapter II for use with general non- 
linear subsystems. In the simulation of a complex nonlinear system, 
if the trapezoidal rule is used for the nonlinear subsystems and the 
bilinear transformation for the linear subsystems, then the same A- 
stable discrete modeling method will have been used for the entire 
simulation. For this reason, the application of the bilinear transfor- 
mation to both rational and irrational transfer functions will be con- 
sidered further in this section. 

If the subsystem transfer function is rational, F^(s) = 
direct substitution of the bilinear transformation s = g(z,h) gives a 
new rational function F'|(z,h) of the same order as F(s). This can 
be used to define a difference equation to be used as the discrete • 
model of the subsystem. It has been observed [3.4] that for high 
order systems, the method is less sensitive to numerical round off 
error if F^(s ) is first expanded into partial fractions and then the 
substitution s = g(z,h) made. A simple example will serve to illus- 
trate the method. 

Example 5.1: Let F(s) be the second order damped system 

F(s) = — 2 & 

s + 2&S+1 
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Substituting s = 


2 ^ z-1 
h * z+1 


there results 


F*(z,h) 



z~l 

z+1 




+ 1 


and expanding 


F*(z,h) = K(h/2) 2 (z 2 + 2z + 1) 

Cl+6h + (h/2) 2 ] z 2 + [-1 + (h/2) 2 ] 2z + [l-6h + (h/2) 2 ] 

If the input to this system is denoted and the output y^, then 
this equation may be interpreted as a difference equation, giving 
the discrete model 


y n = 


[l + 6h + (h/2) 


{K(h/2) 2 Cv n + 2V n _ x + v n _ 2 ] - [-1 + (h/2) 2 


32 y. 


n-1 


- [l-6h + (h/2) 2 ]y n ^ 2 } 

(3.3) 

Note the appearance of v n on the right hand side of equation (3.3), 
which means that the model is externally memoryless, as expected. 

The computed response of this model to a step input is shown in Fig. 

3.1 for different values of the step size, h. For these computations, 

K was taken as unity, and a damping ratio of 6 = .15 was used. The 
error between the exact solution and the computed response for h = .1 
is smaller than the resolution of the graph, hence the curve for h = .1' 
can also be considered as the exact response. 
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If F(s) is not rational, a finite difference scheme cannot be 
used, hence the calculation is done in the time domain, determining 
f(n,h) (the unit response) and using the convolution summation 

n 

y(n) = Y, f(n-k,h)v(k) • 
k=o 

The f(n,h) are identified as coefficients in the expansion 

CO 

F*(z,h) = Y f(n,h)z" n . 
n=o 

Depending on the form of F(s), various techniques can be used to deter 
mine the desired coefficients. One method is a direct expansion in 
powers of z ^ as shown in the following example. 

Example 3.2 : Let 


f(t) = , F(s) = i 

//TTt a/S 

Then applying the bilinear transformation, 



F*(z,h) = / | 

v 1-z v 



♦ ♦ • 


Expanding J 1-z' ^ using the binomial theorem, 



gives 


F*(z,h) = — Cl + z" 1 ' + 1/2 z“ 2 + 1/2 z" 3 + 3/8 z" 4 + 3/8 z 5 + . . . ] 

q/2h 

Thus 

{f(n,h)l =—{1,1, 1/2, 1/2, 3/8 ,3/8, ...j 
«/Zh 


Another method to evaluate f(n,h) involves the use of Laguerre 
functions. For a general subsystem, 


F*(z,h) = F* [g X (s,h),h] = F(s) 


(3.4) 


hence 


F(s) 


I f(n,h) 
n=o 



(1 + sh/2) n 


To carry the development further, a discussion of the properties of the 
Laguerre functions (see e.g. , [3.7], p. 93f) must be given. The 
Laguerre polynomials are defined as 


,n 


l (q) A e 0 *- - — (q n e 
n = dq 11 


= (-1)" (q n - q"- 1 + q n " 2 +... + (-l)V.) 



thus 


V q) 


= 1 


V q) = 1-q 
^ 2 (q) = q 2 - 4q + 2 


The La guerre functions are defined as 




e (q) 

n vn 

nl 


These functions form an orthonormal set, 


I V q) V q) ^ = 6 mn 

o 

and have the generating function 
c P(Q,x) = e 

n=o ' 

A useful recurrence relation for the Laguerre functions is 



(3.5) 


= 


( 2n+l-q) cp n - n c P R _ 1 


n+1 


Their Laplace transform pair is 


(3.6) 


$ (s) = X[cp n (t)] = -^ S - ~ 1/2) n 


n 


(s + 1/2) 


n+1 



Returning now to the problem of finding f(n,h), equation (3.4) is 
expressed in the form 


F(s) = (l + 4r) l f(n,h) 


n=o 


(i - f r 

(x + f r 


or , letting 


F(s) = , and f(t) = 

1 + 


X" 1 CF(s) ] 


as 


f(t) = l f(n,h)(-l) n <p n (f) 


n=o 


Thus f(n,h) equals the product of (-1) and nth coefficient in the expan- 
sion of f(t) in Laguerre functions 9 ^— ) . Note, because of orthogonality, 

00 

f ( n ,h) = (-l) n J f(t) 

o 

Other equivalent relations can be obtained. 

This second method will be demonstrated with an example. 

Example 3.5 : Consider an ideal delay, 


f(t)= u o (t-T) 


F(s) = e 


*sT 
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then 


F*(z,h) 


e 


2T fl-z 




(3.7) 


In this example, effort can be saved if the similarity of (3.7) to 
(3.5) is noted. 

Let 


x = - 




Let 


F(q,x) 





4T\ 

q / 


Thus, from (3.5) 


or 




4T\ 

q / 


1 x n V q) 

n=o 


F*(z, = (1+z 1 ) £ (-D n z' n <f> n (q) 

n=o 

CO 00 

= £ C-D n Z" n cp n (q) + £ (-l) n z" n_1 cp n (q) 
n=o n=o 


(3.8) 
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From (3.8) there results 

f(0,h) = cp o (q) = e"^ 2 

f(l,h) = 9 x (q) -cp 0 (q) = qe ^ 2 

■ 

f (n 5 h) = (-l) n [cp n (q) -tp^q)] n ^ 1 (3.9) 

If the recursion relation (3.6) is used to obtain cp n (q), then f(n,h) 
is easily obtained from the previously computed values of the Laguerre 
functions. It should be noted that the step-size, h, is not limited to 
an integral submultiple of the delay time, T. The model obtained is 
valid for any positive h. 

Plots of f(n,h) for a time delay of three units, T=3, and for 
various values of h are given in Fig. 3.2. These results are difficult 
to interpret, since a discrete model of a unit impulse at T=3 is being 
generated. The plots of Fig. 3.3, however, in which the step response 
of the above model is plotted, demonstrate that the model produces a reas- 
sonable approximation to the true step response, a step at time T=3. 

“S t 

The bilinear transform applied to e results in a more complex 

sh 

model than would the usual z-transform identification z = e . In 
general, consider a transfer function of the form 

F(s) = F(s,e~ sT ) 

where for simplicity, T = kh, k an integer. 
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One might suppose that the transformation 
F*(z,h) = F(g(z,h), z“ k ) 

would have the same stability properties as g(z,h) alone, since 
sh 

z = e also maps the left hand s -plane into the interior of the unit 
circle in the z-plane. This corresponds to using one A-stable trans- 
formation on one part of the problem, and another on a second part. 

““ s t — 1c 

The second transformation, e - z , replaces a delay in the 
continuous time domain by a delay in the discrete-time domain, and 
thus is very easy to implement. The fact is however, that A- 
stability is not generally preserved when different transformations 
are used, as is shown by the following example. 

Example 5.4 : Consider the system shown in Fig. 3.4a, which consists 

of a damped second order system followed by a time delay in a simple 
feedback loop. The second order system has been considered in 
Example 3.2, and the model derived there using the bilinear transform, 
equation (3.3) will be used here. 

If a discrete model for the overall system is formed, combining 
the above model for the second order system, with the simple model 
e = z for the ideal delay, where k = T/h, an integer, the re- 
sultant model will not be A-stable, as can be seen from the following 
argument . 
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(a) CONTINUOUS SYSTEM WITH UNIT DELAY 



(b) MODEL OF CONTINUOUS SYSTEM, USING TWO DIFFERENT 
TRANSFORMATIONS 



(c) CONTINUOUS SYSTEM FOR DETERMINING STABILITY OF 
THE MODEL (b) (ABOVE) 

FIG. 3.4 SYSTEMS FOR EXAMPLE 3.4 
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• — V 

The discrete transfer function z is obtained by usual z- 
_sT 

transformation of e . However, it is possible to consider z to 
be the result of applying the bilinear transformation to the transfer 
function 


G(s) 


, r 1 - sh /2 i 
L 1 + sh/2 J 


G(g(z,h)) = z 


Hence, if the original continuous system were to have G(s) in place 
“ST 

of e as shown in Fig. 3.4c, A-stable modeling (bilinear everywhere) 
would result in the model of Fig. 3.4b. If the system of Fig. 3.4c 
is stable, then the model of Fig. 3.4b will be A-stable. However, it 
is easy to choose values for the parameters for which the original 
system is stable, but for which the systems of Fig. 3.4c and therefore 
the discrete model of Fig. 3.4b are not stable. 

Nyquist plots for both the original system, Fig. 3.4a,- and the 
system of Fig. 3.4c are given in Fig. 3.5, with T=3, K = .39, and 
6 = .15, k=2. The first plot does not enclose the -1 point, but the 
second does. Thus the original system is stable, yet the model of 
Fig. 3.4b, derived using two different modeling techniques for the 
separate parts of the problem, is not stable. 

To verify these conclusions, a computer simulation to obtain the 
step response was performed and the results are given in Fig. 3.6. 

The calculated results for both models using a very small step size, 
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h = .01, coincide. As expected from the Nyquist diagram, the solution 

is a damped sinusoid, settling to some positive value. The result for 

the model of Fig. 3.4b, which mixes two modeling techniques, is shown 

to t be unstable. The remaining curve is the response of the discrete 

model formed by using. the bilinear transformation for both subsystems. 
-sT 

The model for e was derived earlier, equation (3.9). 

It was found that the model of Fig. 3.4b was marginally stable 
at h = .75, and since the natural frequency of the second order system 
is equal to unity, step-sizes larger than .75 are well within the 
limits of the sampling theorem. The model formed using the bilinear 
transformation for both subsystems is, as expected;, A-stable, i.e., 
stable for any value of h. 
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1.0 


-o-o — BILINEAR MODEL, h= 1.5 
-X— x— - MODEL OF FIG. 3.4 b, h = l.5 



FIG. 3.6 COMPUTED 
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CHAPTER IV 


SORTING AND ITERATION 


4 . 1 Introduction 

It was seen in the previous chapters that for both the OSM and 
the ISM methods, a set of equations, usually implicit, must be solved 
at .each time point to provide the simulation outputs. Computational 
procedures for the solution of these equations are presented in this 
Chapter . 
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4.2 System Equations and Directed Graphs 

The OSM method for simulation, Which was presented in Chapter II, 
resulted in a constraint equation (2.6) to be solved at each time 
point t . This equation, in discrete notation, is* 

X n = M GCX^VV f Pu n . (4.1) 

As in equation (2.6), G represents the combined vector of output 
equations for each system described in normal form. The appearance 
of v^ on the right hand side of (4.1) is due solely to direct physical 
connections from input to output in the original continuous system. 

. The ISM method for simulation produced a similar equation to be 
solved, equation (2.12), repeated here for convenience, 


v =MR(v,t)+Pu 
-n v -n’ n -n 


(4.2) 


where R represents the combined vector of discrete models for the 

input-output-state relations of the subsystems, assuming none are 

internally implicit. The appearance of v on the right hand side of 
> *~n 

(4.2) is due either to direct physical connections from input to output 

in the original physical system (which was the only reason that v 

~n 

appears in (4.1)) or to the modeling of dynamic subsystems as externally 
memory less. The case of internally implicit subsystems is left to the 


end of this section. 
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In equation (4.1), the state of each subsystem is assumed to be 
computed before it is needed in the evaluation of each component of 
G, thus x n can be regarded as a known entity. For solution purposes, 
equations (4.1) and (4.2) are thus in the same form, and in what 
follows only equation (4.2) need be considered. 

It will be recalled that the combined matrix [MjP] has only one 
non-zero element, a(+l), in each row, so that equation (4.2) represents 
a sparse system of equations. To attempt solution without taking ad- 
vantage of the structure of the equations would be computationally 
inefficient. For example, consider a simple cascade of N single- 
input single-output subsystems, with an external input, u^, feeding 
the first system. Equation (4.2) for this system is 


— — 


-H 


— 


— — 

V 1 


0 0 0 .... 0 0 




U 1 

CM 

> 


1 0 0 .... 0 0 


wv 


0 

V 3 


0 1 0 .... 0 0 


r 3 (V 3>V 


0 

V 4 

• 


0 0 1 0 0 

. . 1 

I 


wv 

« 

1 

0 • 

• 

_ v n_ 


* • 

0 0 0 .... 1 Oj 


* 

r (v ,t ) 

_ n v n* n _ 

| 

0 

0 


where the time subscript has been deleted on all variables except t 
in order to simplify notation. As one possible approach to solution, 
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consider a simple iteration procedure of the type- used in the trans- 
mission line example of Section 2.4, i.e., 

> 

k+l V 

= M R (v' t) + P u 
— n — n 5 n — n 

which begins with some initial estimate v°. For each iteration, all N 

subsystems are evaluated, and it can be shown that convergence is 

2 

achieved in exactly N iterations. Hence a total of N subsystem 
evaluations must be made. But this would be an extremely inefficient 
method of solution, for the equations can be solved by simply evalua- 
ting them in the proper cascaded order, 


v 

v 


1 

2 


= u. 


= r l (v l> t n ) 


V N = r H-l< v N-l,t > 

J n 

In general, a sequential evaluation of the above type is possible only 
if equation (4.2) is not implicit, i.e., the system is free of directed 
loops. For example, if the input to the first system of the cascade 
connection considered above is connected to the final output, the 
result is a simple feedback loop. If each equation is substituted into 
the succeeding equation, there results, 



86 . 


V N - r N^ r N-l^ r N-2^ * * * * r l^ v N^ * • • ^ ^ 

Thus there exists only one unknown, to be determined. For this 
simple loop, a substitution technique is one obvious method that 
reduces the work to be done in solving the equations. In general, 
however, the structure -is more complex and more sophisticated tech- 
niques will be necessary. 

In order to take advantage of the structure of equation (4.2) 

for general systems, a directed graph representation will be used. 

Each vertex of the graph will represent one subsystem, and the arcs 

will represent components of v^. The external inputs, u^, will be 

considered source nodes. An example is given in Fig. 4.1, in which 

til 

the arc v. . represents the j input to subsystem i. This example 
can be used to suggest a general solution technique based on the 
theory of directed graphs. 

Looking at the graph, it can be seen that there are' no directed 
arcs from the set of nodes [3, 4, 5, 6] to the set [1,2], or equiva- 
lently, there is only feed-forward of information from subsystems 
[1,2] to subsystems [3, 4, 5, 6]. Hence, with v-, = u n (known), one 

X,i j. x 

method of solution is to solve the equations represented by the 
directed loop between systems [l] and [2] and thus to produce v^ 

With v, , determined, the relationships among the subsystems [3, 4, 5, 6] 
can now be solved. 
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Generalizing, it is desired to partition the subsystems (nodes) 
into sequential groups, with each group containing a minimum number 
of subsystems, so that only a feed-forward of information exists from 
each group to the groups that follow. With this partitioning, the 
equations of each group can be solved simultaneously, and the results 
used in solving the equations of the groups which follow. A minimum 
number of systems is desired in each group because this minimizes the 
number of equations which must be solved simultaneously. For each 
group, the subgraph formed by deleting all arcs from the system graph* 
except those interconnecting the nodes of the group will be called the 
group- graph . For complex systems, an algorithm is necessary to accom- 
plish the partitioning and one developed by Steward [4.1], based on 
previous work of Sargent and Westerberg [4.2] is given in Section 4.3. 

The partitioning process is only the first step toward a solution, 

for it now remains to solve the coupled equations of each partitioned 

group. A subset of the variables of the group-graph, called iteration 

variables , are chosen such that it would be possible to evaluate all 

the subsystems of the group in some order if these variables had known 

til 

values. If we denote this subset as w^,- for the i group, then, for 
this group, equation (4.2) can be simplified to 


w. = Q.(w.) 

l l 


(4.3) 


where the function is simply a sequence of subsystem evaluations. 
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If all the subsystems of the group are linear, then equation 
(4.3) is linear, i.e., 


w i = Dw i + d * 

In this ease, the equations of the group can be solved explicitly 
to give. 


w ± = Cl-D] -1 d . 

By definition, the determination of w^ allows the explicit evaluation 
of all the subsystems of the group, and hence the updating process 
for the group is complete . ■ If , however, the group contains any non- 
linear subsystem, then an explicit solution is not generally possible, 
and iteration procedures must be used. 

An iteration procedure is a technique of the form 

w* +1 = EEQ^w^] (4.4) 

k th 

where w^ denotes the k . iterate and the sequence' begins with 
some initial estimate w? . E is some estimation operator which chooses 
a new and hopefully better value for w^ from the previous value. In 
order to reduce the computational effort, the smallest possible set of 
variables w^ should be chosen for each partitioned group. 
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The choice of a minimum number of iteration variables for each 
group can be shown equivalent to a standard problem in the theory of 
directed graphs, the determination of a minimum feedback arc set for 
each group-graph. The following definitions, taken from Younger [4.3], 
will help to make this identification. 

Definition 4.1: For a directed graph, a feedback arc set is a 

set of arcs which, if removed, leaves the resultant graph circuit free. 

Definition 4.2; A feedback arc set is minimum if no other feed- 
back arc set for that graph consists of a smaller number of. arcs. 

Definition 4.3: A sequential ordering of a graph of N nodes is 

a 1-1 function from the nodes of the graph to the integers 1, — ,N. 

Definition 4.4: For a directed graph of N nodes with sequential 
order R, the connection matrix 


C " [C R(i),R(j) ] N,N 

has one row and one column for each node of the graph, and the entry 
at row R(i), column R(j), denoted R(j) = numl)er of arcS from 

node i to node j. 

Younger shows that for any sequential ordering R of the nodes of 
the graph, the removal of those arcs (i,j) for which R(i) ^ R(j) must 
eliminate all directed loops. If all arcs (i,j) for which R(i) s R(j) 
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are selected as iteration variables, then, since no directed loops 
remain, the subsystems can be evaluated in the order determined by R. 
Hence, a sequential ordering determines a choice of iteration varia- 
bles. To find a minimum number of such variables, it is therefore 
desired to find a sequential ordering R^which minimizes the number of 
arcs (i,j) for which R(i) s R(j). This is exactly the problem of de- 
termining a minimum feedback arc set for the group-graph. This 
problem is considered in the next section, but -first a discussion of 
systems containing internally implicit subsystems will be given. 

It was seen in Chapter II that there are two approaches to the 
solution of the equations of internally implicit subsystems. In the 
first approach the equations were to be solved by a process 
internal to the subsystem. Hence, from an input-output viewpoint, 
the equations of the subsystems need never be considered and the de- 
velopment of the proceeding section applies exactly. 

It is possible, however, to treat the equations of an internally 
implicit subsystem simultaneously with the implicit equations of the 
overall simulation. This approach will in general be more efficient 
since the structure of the overall simulation is taken into account. 

A method of accomplishing this will now be given. 

The equations of an internally implicit subsystem can have two 
forms, equations (2.10) and (2.11). Equation (2.11) will be considered 
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first, since the solution technique is .more easily illustrated for 
this form. Equation (2.11) is, ■ 


x n ” q ( x o’* * ’ » x n* v o , ***’ v n’ t o > ** , ’ t n' ) 
y n " r ^ x o’***’ x n ; v o 5 **’’ v n ;t o’* " * 


Since, for this model, the equation is implicit in the state, x , 
an obvious solution technique is to consider the state as an auxiliary 
iteration variable for this subsystem. This in no way affects the 
choice of iteration variables for each group- graph, since only inputs 
and outputs are shown in these graphs. 

If the system equation is of the form (2.10), repeated here. 


y = p(y~ s -.. 5 y > v ,...,v ;t , 
J'n ^ w o’ ,J n’ o > ’ n o’ 


,t ) 
’ n 


then the implicit equation for this subsystem involves the output, y . 
This would represent a self-loop in the group-graph, were it to be added 
there. In order to break all loops, this self loop must be broken. 

Thus y n can be immediately labeled as an iteration variable, and re- 
moved from the group- graph. 
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To remove the arcs corresponding to y^ for a specific 

subsystem, it is necessary only to modify the matrix M as follows. 

til 

The (i,j)th entry of M indicates that the j output is connected 
til 

to the i input. Thus all entries in the column corresponding to the 
output in question need simply be removed. In the following section, 
it will be assumed that this modification has been made, if necessary. 
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4 , 3 Sorting Procedure 

As seen in the previous section, the sorting procedure should act 
in two steps. The first step is to partition the nodes of the system 
graph into sequential groups such that there exists only feed-forward 
of information from any group to those groups which follow it in the 
ordering. This defines which subsystem equations must be solved 
simultaneously and in what order. The second step is to find a mini- 
mum feedback arc set for each group-graph, which determines a set of 
iteration variables for each group of equations. 

A preliminary simplification of the system graph is possible 
which reduces the complexity of the partitioning process. The source 
nodes (external inputs) and the arcs which leave them, can never be 
iteration variables, for they represent known quantities. In es- 
sence, the source nodes can be considered as a first group of the 
partition, for their "equations” are always solved by definition. 

Hence in what follows, the matrix of external input connections P," 
need -not be considered. 

In order to apply the partitioning algorithm, the connection 
matrix, C, must be determined. All the information necessary to 
determine this matrix is contained in the matrix M of equation (4.2) 
from which the system graph was derived. Adding together all the ' 
rows of M corresponding to subsystem i, i=l,...,N and adding together 
all the columns corresponding to subsystem j, j=l,...,N, there results 
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a new matrix, S, whose (i,j) element is the number of arcs from 
subsystem j to subsystem i. Hence, the connection matrix is the 
transpose of S, or C = S t . 

Once the connection matrix is determined, the algorithm presented 
by Steward [4.1] can be used to partition the subsystems into the 
desired groups. This algorithm proceeds as follows: 

(1) A column is sought having no predecessors, i.e., no entry 
in the column, and that column and the row corresponding to it elimi- 
nated. This process is repeated until there are no further columns 
without predecessors. 

(2) Any path traced following predecessors must lead into a loop. 

Since the graph is finite and by Step (1) there is no column without 

a predecessor, such a path must eventually repeat a column. Retracing 
the path from this column gives a loop. 

The union of two columns is defined to. be the sum of the columns, mod 2. 

(3) When a loop is found, the set of columns in the loop is 
replaced by one column which is the union of the columns replaced. 

This is called collapsing the columns of the loop. Similarly, the 
rows corresponding to these columns are collapsed. Any entries which 
result on the diagonal are set equal to zero, since they simply re- 
present connections within the collapsed loop. . 
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When a column without predecessors is eliminated, that column 
and the columns which collapsed to form it represent the subsystems 
in a group. The order in which columns without predecessors are 
eliminated gives an order in which these groups may be solved. 


Example 4.1: The connection matrix for the graph of Fig. 4.1, 

eliminating the known input v^ 1 and source u^, is 


C = 


1 

2 

3 

4 

5 

6 


1 

0 

1 

0 


2 

1 


0 0 
0 0 


3 

0 


1 
0 
0 
1 

0 0 


4 

0 

0 

1 


1 

0 


5 6 

0 o' 


0 

0 


0 0 


0 

1 

1 


0 0 


0 


There are no columns without predecessors. Beginning in column (1), 
the loop 1-2 is found. Collapsing the loop gives 


1-2 

3 

4 

5 

6 


1-2 3 

'0 1 
0 0 
0 0 
0 1 
0 0 


4 

0 

1 

0 

1 

0 


5 

0 

0 

0 


6 

o" 

1 

1 


0 0 
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Column (1) now has no predecessors, and .can be removed. Thus, [1,2] 
is the first group. Removing this column and row gives 

3 4 5 6 

3 To 1 0 1 

4 0 0 0 1 

5 110 0 

6 1_0 0 1 0 

The loop 3-5-6 is now found, and collapsing these rows and columns 
gives , 

3-5-6 4 

3-5-6 p 0 1 

4 1_ 0 

These two columns form a loop, hence the second group is [3,4, 5, 6]. 

This grouping was previously obtained intuitively. 

With the partitioning accomplished, a minimum set of iteration 
variables must be found for each group. In the previous section, it 
was stated that the problem of determining a minimum set of iteration 
variables corresponds to the problem of determining a minimum feedback 
arc set for the group- graph. In order for this correspondence to hold 
for all systems,- a simple modification must be made in the group-graph, 
under certain conditions, to adjust for the following difficulty. In 

/ 

the group-graph, a single output of one subsystem going to many inputs 
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represents only one variable, but many arcs. However, if the subsystem 
has separate outputs going to each of the other inputs, then this 
represents - many variables and as many arcs. The two graphs look 
identical except for labeling of the arcs, but are very different as 
far as choice of iteration variables is concerned. If the single 
variable of the first situation is chosen for iteration, then this 
removes all the arcs which it labels. To break all the arcs in the 
second case, all the variables must be chosen for iteration. To 
remedy this difficulty, a single output going to a multiplicity of 
inputs will first be removed from the node through a single arc to a 
new inserted node which' will then branch to the necessary input loca- 
tions. The cutting of the single arc connecting the original node and 
the inserted node will always be equivalent to the cutting of all the 
branching arcs. 

Once the graph has been modified as above, a minimum feedback arc 

set for the group-graph is a minimum set of iteration variables. It 

\ 

was seen in the previous section that an ordering of the nodes deter- 
mines a feedback arc set. Hence, such an ordering will be determined 
according to the following definition [4.3]. 

Definition 4.5 : An optimum ordering R is a sequential ordering 

of the nodes of a directed graph for which {(i,j), R(i) s R(j)} is a 
minimum feedback arc set. 
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Before an optimum ordering is sought, recognition of certain charac- 
teristics of the graph may simplify the problem. The first concerns loops 
of two arcs. A theorem of Younger’s gives this simplification. 

Theorem 4.5: The set of optimum orderings for a given graph is 

invariant under the removal of loops involving two arcs. In the 
case of multiple arcs between nodes, e.g. , p arcs from node i to node j 
and q arcs from node j to i, with p > q, then up to q arcs may be re- 
moved from i to j and the same number removed from j to i. 

Hence, when searching for an optimum sequential ordering, all loops 
involving only two nodes can .be removed from the graph. As an example, 
the first group of the graph in Fig. [4.1], nodes [1,2], contains a 
loop which can be removed. In this case, both possible orderings are 
optimal, for the ordering [1,2] gives v, 9 as the only iteration vari- 
able, and the ordering [2,1] gives as the only iteration varia- 

ble. 

A second simplification, which is an extension of a simplification 
first given by Sargent and Westerberg [4.2], can be accomplished if the 
form shown in Fig. 4.2a appears in the graph. The arcs to the left 
and right of nodes 1 and 3 respectively denote connections to the 
remaining nodes of the graph. For this configuration, to cut all 
paths 1-2-3, only' the min[p,q] arcs need be cut. Node 2 therefore 
can be temporarily removed, and the arcs from nodes 1-3 replaced, 
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b ) 

POSSIBLE ORDERING 

1 2 3 

2 1.3 

1 3 2 

3 2 I 

2 3 I 

3 I 2 


NUMBER OF FEEDBACK ARCS 


s 

p + s 
q + s 


NODE I 
► PRECEDES 
NODE 3 


p + q + r 
P + r 

q + r 


NODE 3 
► PRECEDES 
NODE ! 


FIG. 4.2 SPECIAL CONFIGURATION WHICH LEADS TO SIMPLIFICATION 
OF THE SYSTEM GRAPH 
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using Theorem 4.5, by a number of arcs r-s+ min[p,q], where a positive 
result is interpreted as arcs from 1-3 and a negative result as arcs 
from 3-1. If the chart in Fig. 4.2b is examined, it is seen that the 
final position of node 2 is determined by the relative position of 
nodes 1 and 3. If 1 precedes 3, then the final ordering should be 
1-2-3 and if 3 precedes 1, then the final ordering should be 2-3-1 
if p < q or 3-1-2 if q < p. 

The simplification techniques presented above greatly reduce the 
effort involved in finding an optimum ordering; however, the algorithm 
developed by Younger will work whether or not they have been applied. 
Younger’s algorithm is given in complete detail in [4.3], but the 
* salient characteristics are presented here. 


Definition 4.6: For a sequential ordering R of a directed 

graph, a consecutive subgraph is any nonempty subgraph composed of 
nodes consecutively ordered by R and the arcs connecting them in the 
overall graph. 

For subgraphs and of a directed graph G which contain 

no common node, let C„ „ denote the number of arcs each of which 

G 1 G 2 

leaves a node in G^ and enters a node in G 


Definition 4.7: An ordering R for a directed graph is admissible if 

(a) the graph ordered by R satisfies C„ _ ^ C„ „ for all appro- 

G 1 G 2 G 2 G 1 

priate consecutive subgraphs G-^Gg 
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(b) the feedback arc set determined by R contains no proper sub- 
set that is also a feedback arc set. 

An admissible ordering for a group-graph is intuitively a good 
ordering. A first admissible ordering is easy to find and in fact, 
the search process could be stopped here if it were decided that a 
small, but not necessarily minimum, number of iteration variables 
would be acceptable. Younger’s algorithm begins by finding an ad- 
missible reference ordering. A new ordering R T , which has a smaller 
number of feedback arcs, is sought, through a perturbation process. 

This new ordering R f becomes the new reference ordering, and this 
process repeated until a reference ordering is obtained which cannot 
be improved. This is an optimum ordering, and determines a minimum 
number of iteration variables. 

Example 4.2: A brief sketch of the application of this algorithm 
to the graph of Fig. 4.1 will now be given. For the second partitioned 
group, nodes [3,4,5,6], the first step is to choose an initial ordering 
and check to see if it is admissible. If it is not, the testing 
process yields the information needed to change to an ordering which 
is admissible. For the example, the given numbering of the nodes, 

[3, 4, 5, 6] serves as a first try. This is shown to be not admissible, 
because the subgroup [3,4] has no arcs going to the group [5], while 
[5] has two returning to [3,4]. The algorithm permutes the nodes to 
obtain the admissible reference ordering [3, 6, 5, 4]. This ordering 
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gives the two iteration variables v^ ^ and v^ ^ • It is seen in the 
graph that if these two arcs are removed, no loops remain. However, 
this ordering is not optimum, for the algorithm next produces the 
ordering [5,3,4,63 with only one feedback arc, v,- It is seen in 
the graph that removing only V(- , does break all loops. Fig. 4.3 

D,J- 

shows how the overall . simulation would proceed, once the sorting in- 
formation is available for this example. 

Younger’s algorithm has one characteristic, already alluded to, 
that makes it more suitable to computer sorting of the equations than 
other techniques. This algorithm can be stopped at any point, and 
the results used meaningfully, because whatever admissible refer- 
ence ordering has been found at that point serves to determine a good 
(but not optimum) set of iteration variables. The other applicable 
techniques, in particular dynamic programming, must work completely 
to a solution before any choice of iteration variables is possible. 
For a large problem, .it may be preferable to preset an amount of com- 
putation time that will be allowed for this stage of the simulation, 
and stop the search if this time is reached. With Younger’s algori- 
thm, the simulation proper could be begun at this point, whereas with 
dynamic programming the search has simply failed. 
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FIG. 4.3 FLOW DIAGRAM OF SIMULATION FOR EXAMPLE 4.2 
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A summary of the entire sorting procedure is given below: 

(1) Partition the nodes with the partitioning algorithm. 

(2) Modify any outputs which connect to more than one input to 
include an extra node. 

(3) Apply the two simplification methods to each group-graph 
(optional) 

(4) Find a sequential ordering for each group which determines 
a minimum set of iteration variables (optionally, a suboptimal set). 

The final result is the sequential ordering of the groups, with each 
group properly ordered (and thus determining the iteration variables). 



106 . 


4.4 Iteration Procedures 

In the preceding section, the equations of each partitioned group 
were shown to have the form of equation (4.3), repeated here 

w i = ^i^ w i^ (4.5) 

where Q_^ is a sequence of subsystem evaluations determined by the sequen- 

til 

tial ordering, for the i group-graph. The ordering was found 
which minimized the number of variables, w.. In this section, techniques 
for solving (4.5) for each group-graph will be discussed, where it is 
assumed that the group contains at least one nonlinear subsystem. The 
purpose of this discussion is to show how the structure of the system 
model influences the choice 'of an iteration procedure. 

In order to simplify notation, the subscript i will be dropped in 
what follows, since only the equations of one group will be considered 
at a time. It will be recalled that equation (4.5) must be solved at 
each time point, and that the time subscript, n, was previously removed. 
Thus in what follows, the solution of equations of the form- 

w = Q(w) 

are considered. 

The first method which suggests itself is direct iteration of the form 

k+1 k. O rn a\ 

w. = Q(w ) , w (4.6) 


Throughout this section the superscript notation ( 
"k-th iterate" . 


indicates 
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This method was successfully used in the transmission line example 
of Section 2.4. The properties of this iteration procedure follow 
easily from some definitions and theorems of analysis (see, e.g. , 
[4.4], p. 627ff ) . 

Definition 4.8: A mapping Q(w) of the w-space into itself is 

said to be a contraction if there exists a c < 1 such that 


j|Q(w) - Q.( w T ) |1 = c|jw-w l jj for all w,w f 


Here H»|| denotes a vector norm. 


e.g. 


ii a max 
x|! 4 i 




if 


Theorem 4.6: Given a continuously differentiable function Q(w), 


||Qj| = c < 1 for all w, 


then Q(w) is a contraction. Here ||.|| denotes the matrix norm 


ilo.ll a m f s|q,J- 

i 
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Theorem 4.7 : If Q is a contraction mapping, a unique solution 

w* of equation (4.6) exists, w* can be obtained as the limit of a 
sequence where 

w k+1 = Q(w k ) k = 0,1,2... (4.7) 

and w° is any given initial value. 

k 

The rate of convergence of w to the solution is given by 



The direct iteration procedure thus converges only if the norm of 
the Jacobian of Q satisfies a Lipschitz condition. The Lipschitz con- 
stant, c, determines the rate of convergence of the sequence and if 
c > 1, the sequence may diverge. For most systems, the function Q is 
not a contraction and more sophisticated iteration procedures must be 
used. 

At this point it is worth noting the effect of integration step- 
size on the convergence of the iteration procedure. The implicit equa- 
tions will often result from the use of A-stable modeling of subsystems, 
in which case the step-size h will appear as a parameter in the func- 
tion Q. In this case, Q. will often be a contraction only for values 
of the step-size less than some value, h c . If direct iteration is to be 
used to solve the equations, then to achieve convergence, it must be re- 
quired that h be less than h c .' But this reintroduces a limitation on h. 
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which negates the purpose' of A- stable modeling. Thus direct iteration 
is not generally useful when A-stable modeling has been used. 

Another useful technique is Newton-Raphson, which for equation 
(4.6) is 


w 


k+1 


= w k - [i - pjV [w k - Q(w k )] k = 0,1,2, , 


where 


^ ^ l w _ k 


The difficulty in applying Newton-Raphson to the solution of 

equation (4.6) for general systems results from the necessity of 

k k 

evaluating the Jacobian at each iteration. In general, may be 

very difficult to compute, or may simply not be available. If the 

Jacobians of each individual subsystem are included in the mathematical 

model of that subsystem, then the sorting procedure of the previous 

section can be used to select the proper derivatives as the loops are 

traversed. It would be desirable however, to avoid as much of the 

aforementioned computation as possible. A method due to Broyden [4.5] 

has this property. The method of Broyden chooses a new estimated value 


for 


[>-<] 


-1 


from the values of the functions at each iteration. The 
method uses this estimated value in place of the true value in a Newton 
iteration, and hence belongs to a class of methods called quasi-Newton 
methods . 

This method is outlined below, with slight modifications due to 


the special nature of the simulation equations being solved. 



no. 


(1) Obtain an initial estimate w° of the solution. The final 
value for w at the previous time point will usually be a good choice. 

(2) Obtain an initial estimate of the iteration matrix 

m ' 1 = h °- 

Again, the final value for H from the previous time point is a very 
good choice. 

(3) Compute f^ = w k - Q(w^) 

(4) Compute p^ = -H^f^ 

k 

(5) Choose t such that if 


then 


w 


k+1 


= w 


k.k 
P t 


Hf k+1 n < iif k n . 


Broyden shows C4.6], that if the method is close to a solution of (4.5), 
then the Newtonian value of unity for t can be chosen and 
the norm minimization avoided. Since, in the problem considered here, 
good initial estimates are available at each time step, t will be set 
to unity. 

k+1 

(6) Test f for convergence. If convergence has not been reached, 

(7) Compute y k = f k+1 - f k 

(8) Compute I? +1 = H*- - (H k y k - p k t k )(p k ) t H k /(p k ) W 

(9) Repeat from step (4). 
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Broyden's method avoids the computation of the Jacobian, but the 
cost may be an increased number of iterations at each step. However, 

t 

since a good initial estimate of the solution and iteration matrix are 
available at each step, in general the number of iterations will not be 
much higher than that achieved using true Newton -Raphs on. 

The following example illustrates the necessity of using the 
more complex iteration procedures. 

Example 4.2: Consider the transmission line example of Section 

2.4, with the capacitor C 2 removed from the nonlinear termination 
(see Fig. 2.2). With the capacitor removed, there results a term in 
the Jacobian corresponding to the incremental resistance of the 
right hand termination, i.e.. 



In the linear case, this term equals R^, and with a value of 2K, the 
norm of exceeds unity in magnitude and it is expected that direct 
iteration would diverge. Similarly, for the diode load, 


9v„ 

^ 

3i, 


V T 


which for small values of i 2 also exceeds unity in magnitude. It was 
confirmed experimentally that direct iteration diverged for either the 


linear or diode load. 
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Both Newton-Raphson and Broyden* s method were now applied to the 
problem with good results. For the linear load, with h = 1/32 and 
€ = .001, both methods' took almost exactly the same average number of 
iterations per step, 1.78. This is not surprising, for once Broyden* s 
method has solved the equations at the first time point, it has con- 
verged to the exact value of the iteration matrix, which is constant in 
time for the linear loads. Hence, after the first time point, there 
is essentially no difference between Broyden* s method and Newton-Raphson. 
In the nonlinear case, both methods converged, with Broyden’s method 
taking a slightly higher average number of iterations per step, 1.92 
for Broyden’s versus 1;83 for Newton-Raphson. Computationally however, 
this means that Broyden* s method was more efficient, since two extra 
function evaluations are needed to compute the necessary derivatives 
for Newton-Raphson. 

The above analysis and computational results suggest that it is 
useful to have available a good hierarchy of iteration schemes, using 
direct iteration if possible, followed by Broyden T s method if the 
former diverged. Ideally, however, a general purpose simulation 
language would have a large repertoire of iteration routines, similar 
to the presently available repertoires of integration methods, with 
the choice left to the user. 
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CHAPTER V 
CONCLUSIONS 


The primary goal of this research has been to develop computational 
methods for the digital simulation of continuous systems. It was de- 
sired to develop methods which would be applicable to a wider class of 
systems and would be more efficient than currently available simulation 
languages . 

By analyzing current simulation methods, it was possible to dis- 
tinguish and describe two fundamentally distinct approaches to simula- 
tion, herein named the OSM and the ISM methods. The OSM method, in 

' # 

general use in available simulation languages, models the overall sys- 
tem, its component parts or subsystems having already been incorporated. 
The ISM method models the subsystems and then interconnects these 
models. It was shown that the ISM method is applicable to a wider 
class of systems than the OSM method. Thus, to achieve greater 
generality for a simulation languages, the ISM method should be used. 
Within this context the types of numerical modeling techniques admis- 
sible for the ISM method were determined, and the characteristics of 
the simulation. which result from using these methods were investigated. 

One problem which limits the generality of available simulation 
languages is their inability to solve sets of implicit equations. 
Previously, this problem was not considered a limitation serious enough 
to justify structuring a simulation language around the solution of 
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implicit equations. As. was shown in this research, however, the solution 
of implicit equations is related to another very serious limitation of 
current simulation languages: their inability to simulate stiff systems 
efficiently. 

A-stable modeling is discussed as a satisfactory technique for the 
efficient simulation of stiff systems, and a new development for the 
application of A-stable methods to modeling of linear systems is given. 
It was shown that A-stable modeling results in implicit equations for 
the overall simulation. Thus, since the practical, efficient handling 
of stiff systems involves implicit equations, solving these equations 
becomes a necessity and justifies the structuring of simulation lan- 
guages around this problem. 

In order to structure simulation languages to deal with implicit 
equations, sorting and iteration procedures for the efficient solution 
of these equations are presented. Subsystems are sorted to obtain a 
good computational sequence. Since nonlinear systems are generally 
present, iteration procedures to solve the system equations are 
presented. 



